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FOREWORD 


This document presents the Final Report on the development and use of 
a computer program for "Hydroelastic Modal Studies" for George C. Marshall 
Space Flight Center (MSFC) of the National Aeronautics and Space Administration. 

This capability was developed as program additions to the NASTRAM (NASA 
Stru cture ^alysis) computer program and delivered to MSFC tor use in the 
analysis of the Shuttle ET tank hydroelastic modes. 

Messrs. D. N. Herting, R. L. Hoesly, and D. L. Herendeen of Universal 
Analytics, Inc. were the primary contributors to the project. Mr. R. McComas 
of MSFC monitored the project and also contributed to the testing of the system. 
Messrs. D. Harper and J. Moreman of Computer Science Corp. performed valuable 
aid in implementing the system at MSFC. 

The Final Report is divided into three volumes. This, the first volume, 
describes the theoretical basis of the system and results of the test cases. 

Also Included are a brief summary of the project history and explanations 
and conclusions based on the results obtained in the contract. The second 
volume contains the Program Documentation and the third volume serves as a 
guidebook CO the use of the system. 
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1.0 INTRODUCTION 


This report describes the technical effort performed in the implementa- 
tion of a genertd three-dimensional hydroelastic capability in the NASi'RAN 
(MSA ^uctural Malysls) computer program. Although NASTRAN had provided 
capabilities for the analysis of compressible fluids with axlsymmetrlc 
geometry, the Space Shuttle External Tank Analysis required more general 
capability and more efficient solution procedures. The basic approach 
described in this report extends the capabilities to provide for arbitrary 
fluid shapes, including tilted free surfaces, and allow for more efficient 
methods of obtaining the solutions. 

One goal of the program development was to provide a general method 
for analyzing the combined mode saapes of arbitrary fluid and structure 
finite element models. The fluid is modeled with three-dimensional solid 
elements with options for tetrahedron, wedge, and hexahedron shapes. The 
elements are connected to fluid grid points which define the pressure in the 
fluid at the specified location. The structure may be modeled arbitrarily 
using the existing NASTRAN elements. The fluid/structure Interface and the 
free surface are defined by the user with special NASTRAN boundary elements. 
A special purpose mesh generator program has been provided to generate the 
actual NASTRAN data cards for the fluid, the structure, and the boundary 
elements for typical tank-type models. 

A second goal of the project was to formulate the matrices to provide 
efficient solutions for large-order problems. The structure matrices are 
processed separately and may be reduced using matrix condensation procedures 
(0MIT) or through a modal formulation using the normal modes of the empty 
structure as solution coordinates. The fluid matrices are then transformed 
and connected to the reduced structure coordinates resulting in small, sym- 
metric, solution matrices. This approach is particularly valuable when 
several different fluid levels are to be analyzed for one structure. The 
structure formulation and reduction Is processed only once. The additional 
calculations for each different fluid case require only fluid matrix opera- 
tions and solution processing. 
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User convenience was provided in the system with the implementation of 
several alternate solution paths and modeling options. These options » which 
allow a wide variety of problem types and provide the user with efficiency 
and accuracy trade-offs, are summarized below. 

• The Direct formulation option uses structure grid point coordinates 
and free surface displacements as solution degrees of freedom. The 
structural matrices may be reduced using the NASTRAN matrix conden- 
sation technique (ASET or 0MIT data) for more economical process5.ng 
of large-order problems. 

s The Modal formulation option uses the mode shapes of the empty 
structure as generalized solution coordinates. With this option, 
the structure may be represented accurately by relatively small 
matrices. 

• Gravity effects are provided which affect both the free surface 
displacements and the structure-fluid Interface. Free body motion 
of the whole system is not constrained as in seme other methods. 

The gravity effects may also be Ignored by simply providing 
single-point constraints on the free surface pressure coordinates. 

s Symmetric systems with symmetric or antl-symmetrlc solution cases 
may be modeled by providing single-point constraints on structure 
displacement and fluid pressure degrees of freedom. 

• Compressibility effects may be provided for cases in which the fluid 
is completely enclosed. (Any constraints on the fluid pressure 
coordinates will act as an opening.) A factor is provided to define 
the overall pressure versus volume change. An alternate is provided 
to constrain the volume change to zero for pure incompressibility. 

• Restart logic is provided in the DMAP (Direct Matrix Abstraction 
Program) to allow changes in the fluid model without reformulating 
the structure matrices, or generating structure modes. A user 
parameter provides this control, Independent from the NASTRAN logic. 

All of the above capabilities were specifically designed for the large- 
order finite element models anticipated for use in the analysis of the Space 
Shuttle tanks. In many cases, the actual computer hardware availability 


would be a significant factor in solving the real structures. Thus, the 
design contained the alternate paths to obtain either accurate results with 
long run times or less accurate results with shorter run times and faster 
turnaround* 

The effort for the project was subdivided into five major tasks, each 
defined by delivery items and milestones. Brief descriptions of the tasks 
and the problems encountered are summarized below: 

TASK I - Preliminary Specifications : The complete theoretical design, 

user Inputs, and program flow was defined in this documentation. 

The report was delivered on schedule within one month of the con- 
tract start date of March 1, 1976. 

TASK II - Phase 1 Program : The basic hydroelastic capabilities were 

installed in a Levell5.9 version of NASTRAN. The program was 
installed at MSFC on the IBM 360/55 computer on schedule during the 
week of April 19. The basic capabilities Included the fluid elements, 
free surface elements, and non-overlapping fluid-structure boundary 
definitions. Results were obtained for the hemispherical tank test 
problem. 

TASK III - Phase 2 Program : The final version of the program was 

installed on the Level 16.0 version of NASTRAN. The delivery, on 
August 9, 1976, to MSFC was delayed approximately one month due to 
delays in the delivery to MSFC of this version of NASTRAN by the NASA 
distribution center. 

TASK IV - Demonstration and Training : This task consisted of three 

parts - (1) demonstration of the system using small problems with 
known solutions, (2) instructing MSFC personnel in the use and 
maintenance of the system, and (3) assisting MSFC personnel in the 
modeling and execution of the large-order Space Shuttle tank problems. 
The first two tasks were performed at MSFC during the week of 
August 9, in conjunction with the program delivery. The third task 
continued through the remainder of the contract. 

TASK V - UNIVAC 1108 Conversion and Final Report : The original schedule 

assumed that the 1108 version of NASTRAN Level 16 would be available 
:o NASA by October 1976. The IBM hydroelastic program would then be 
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converted to the 1108 at MSFC by UAI. However, the actual receipt 
of the 1108 Level 16 NASTRAN did not occur until January 1977, 
delaying the project by three months. The final report, wnich 
consists of program documentation, awaited the 1108 conversion and 
the execution by MSFC of the large-order Shuttle tank problems. 
However, preliminary versions of the user documentation were 
delivered to MSFC in August 1976, thereby allowing knowledgeable 
use of the program. This final report, contained herein, completes 
the outstanding contract deliverables. 

The theoretical approach and the detailed analytic steps used in the 
formulation of the problem are described in the following theoretical 
approach chapter. The results of the test and demonstration effort are also 
summarized in this volume of the report. Volume II of this report contains 
the descriptions of the user-supplied data to the program and instructions 
for its use. Volume III contaitis the programmer's reference material 
describing the new code and modifications to NASTRAN. 
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2 . 0 THEORETICAL DEVELOPMEMT 


In this section the theory Is developed for the ,,enerel thtee-dinenslonal 
hydroelestlc analysis In kasIRAN. Both the structure and the Interacting 
fluid will be Ideallaed as general, three-dinenslonal Unite element models. 
Effects of f.ee surfaces and steady-state gravity will be included. Th» 
fluids are assumed to be incompressible, irrotatloual. and non-viscous. 

Small motions of both structure and fluid .eiac.'ve to the static solution • 
will be analyzed. 

The basic development of the finite .’.meat creations for small motions 
of fluids is described in Reference 1. In Rcr^.ence 2. the basic equations 
are cast in the form of integrals representing the time derivatives of Energy 
and Work using the fluid pressures as the unknown coordinates. The scalar 
pressures, rather than three displacements, will be used as degrees of freedom 
at each point in the fluid, which avoids extraneous rotational motions and 
directly provides for Incompressibility. The disadvantage is that the struc- 
ture and fluid are not automatically connected at their boundary. The 

pressures in the fluid must be related to the displacements of the boundaries 
through area factors and flow relationships. 


2.1 FLUID FIELD EQUATIONS 

In Reference 2 the fluid field equations ere developed lu the for. of 
energy iutegrele using principles of variational calculus. The basic result 
for the compressible case is the equation: 


®[/v (4 P • ’p) ■•v] - I ip(j 7pj ■ dS . 0 (1, 


where : p 

P 

3 

P 

V 
dS 

V 
6 


is the pressure 

is the time derivative of pressure 
is the bulk modulus 
is tha mass density 
is the volume 

is an incremental surface vector (outward) 
is the vector gradient operator 
is the variational operator 
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With the incompressib?^e fluids, the bulk modulus is assumed infinite and 

.2 

the p term disappears. On the exterior surface the pressure gradient may be 
replaced with the acceleration vector using the basic momentum equation: 

Vp = - pii (2) 

Equation (1) therefore becomes: 

6 f (^P • ^P) dV + I 6pu ♦ dS = 0 (3) 

or 

6U + 6W = 0 (4) 

In the finite element method of solution, a set of variables, equal to 
the value of p at specific points, is chosen and the volume Is divided into 
subregions, called fluid elements, with vertices defined by the location of 
the variables. As with structural finite elements, a "stiffness” matrix [K] 

Is formed frcm the potential energy U by the equation: 


ij 3Pi 8Pj 

The "equivalent" potential energy for each subregion is, from Eq. (5), 


(5) 


u = f ^ Vp . Vp dV 
J V ^ 


( 6 ) 


The pressure field for each subregion (fluid element) Is dependent on the 
pressures p^^ at its vertices. 

The second term of Eq. (5) defines the "load" vector, {Q}, on the boundary, 
which excites the fluid. For each point j on the boundary the load'ng value is 


Ql 


3(6W) 

3(^ 


(7) 


Since the accelerations, U, are defined by the structure, a transformation 
matrix may be developed from Eq. (?) such that: 


a 
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^Q> - [BHii} 


where 


( 8 ) 




3p. 3u, J P“ dS 

J 


(9) 


Note that the units of B 


squared. 


are area and the unite of Q are volume 


per time 


If we combine Eqs. (4) (fk\ rj\ j ,ov 
relating the structure and’fluil cco:dinates 


IK ]{p) + [B]{a> = 0 


On the other hand, the fluid effects the 

the structure surface area The 1 ^ ®PPlyri»g forces over 

area. The incremental w^rk done on the structure 


( 10 ) 


Is: 


6W 




P(6u • d') 


( 11 ) 


Iha force on each point atructure degree of freedom, F ia obtained f 
work by the equation ’ j * obtained from the 


equation: 


3(6W ) 
F. = ® 


j 3 (6uj ) 


( 12 ) 


PotnT.To, diaplacementa „ are function, of the grid 

.oln dU lerement, u Equations ( 11 , and ( 12 ) are evaluated resulting in 
the transformation matrix. A. where: resulting in 


= lA]{p} 


(13) 


and 


ij 


du 




pu dS 


(U) 


Comparing Eqs. (9) and (14) „e observe the 


L\t 
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[A] = [B^] 


s s 

If [M ] and [K ] are the mass and stiffness matrices of the structure, the 
matrix equation for the structure coordinates is 


or, from Eq. (12): 


+ [K®]{u} - {F} - {0} 


+ [K®Hu} - [A]{p> = {0} 


Equations (10) and (17) become the system of equations for a solution. 

2.2 FINITE FLUID ELEMENTS 

Three types of fluid elements are used to represent the three-dimensional 
fluid; the 4-point tetrahedron, the 6-point "wedge," and the 8-point hexahedron, 
llie wedge and hexahedron elements are fabricated from three and ten tetrahedra, 
respectively, as ihown in Figure 1. The pressure function within each tetra- 
hedron is assumed to be a linear function in three directions, or; 

P = q^ + q^^x + q^y + q,z (18) 


The transformation between the pressures at the grid poincs P^, i = 1,2,3 
and 4, and the coefficients q may be obtained from the matrix equation: 


{p} = 


/% 


1 *2 >'2 '2 
1 Xj yj z. 

1 ^ h 


The matrix in Eq. (19) may be inverted, giving the matrix H, where 

{q} = [H]{p} 


"k ■ 5»KjPj 








STRUCTURAL ELEMENTS 




(a) Tetrahedron. 




(b) Wedge and One of Its 3ix Decompositions. 



(c) Hexahedron and its Tv/o Decompositions. 


FIGURE 1. THREE-DIMENSIONAL FLUID ELEMENTS 


The vector gradient of the pressure is obtained from Eq. (18), giving: 


Vp = qj^i + q^j + q^it 


( 22 ) 


From Eq. (5) in the previous section, the energy function, U in the element 
is: 




From Eq. (6), the stiffness matrix terms for the element are: 


(23) 


V . 


2 

3^U 




Using the "cha?.n" rule for differentiation, we obtain 

K - "S' 

M 3^2 SPi SPj 

and from Eq. (21) : 



k=l 


p ^k)l\i^S,i 


(24) 


(25) 


(26) 


(Note : 0 if k # A and « 1 if k " Z) 

Equation (26) may be cast as a matrix product defining the fluid "stiffness" 
matrix for the tetrahedron, IK^], as: 


[K^] - ^[h'^][H] (27) 

The wedge and hexahedron elements are generated by adding the appropriate 
stiffnesses for the component tetrahedra. 

2.3 FLUID/ STRUCTURE BOUNDARY MATRICES 

As defined in Eq. (14) of the general development, the area matrix [A] is 
defined as 
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pu dS 


(2 


3u^ 9p^ 

where u and p are the displacements and pressures at the surface, S. The 
intersecting areas of the structure and the fluid are specified by the user 
as fluid-structure element pairs. From elementary geometry, the locations of 
the fluid points and the structure points are obtained in a coordinate system 
on the fluid face. Equation (28) is evaluated for each intersecting area of 
structure and fluid. For simplicity, only triangular structure elements are 
considered below. Quadrilateral elements are treated as four overlapping 
triangles . 

Several possible examples of overlapping areas are shown in the sketch 
below. 



Clearly, the number of combinations is too numerous to identify each case 
and provide a specific set of equations. Rather, a general algorithm will be 

developed below. The basic logical steps, listed below, will apply to all 
cases* 



thTn!iV:::r" 

points e whIreT”^U^^^ outside the area, search for any intersection 
Foxncs, e^, where two lines cross. 

r:rr -- — - 

.he ej : :::: ziZ"7 1 

a fatal error has occurred. 

If inrersecrloee occur, the overlapping area la detemlned aa follca. 
Starting with a cotnaon ed^a nn-fn#- *.u •* 

rr.:- r t- 

Iculated. An eaanple la aho«n in the aketch below. 



The .lat of polnta.a. defining the area ares 

“ structure point 
6 j^ - edge point 
^2 ^ point 

^2 ~ fluid point 
®2 “ ®*^ 8 e point 

The area of the polyhedron la obtained fron the line Integral: 



( 29 ) 
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■ tfiii iiiiiri 


JiSm 




y gw'i' wm, ly lui lu ) I j 


5. Using the finite element displacement functions, the coefficients C .(x ,y ) 

are evaluated at all points on the polyhedron. The pressure applied^o^ach 
point in the area is : 


P 


s 


2 C ,p. 


(30) 


where is the pressure on each polygon point and p. is the pressure at 
corner j of the fluid element. 

fluid elements the pressure distribution is 
P = + q^x + q^y (31) 

Evaluating the equation at the four corners results in the formula: 



The coefficients C(x,y) are therefore: 


(32) 



(33) 


On quadrilateral areas, an isoparametric distribution of pressure is used. 
The isoparametric coordinates are shown in the sketch below. 
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i 

4 


3 



( 


The pressure distribution is; 

" ■ Vl + f2P2 + V3 + V4 

where; . (1 - 5)(1 - n) 

^2 = ^(1 - n) 

= a - On 


and the corresponding location definitions are 


x(C,n) 

y(C,n) 


*1*1 ^2*2 V 3 + ^4*4 

Vl + V 2 V 3 <^4>'4 


If and for a point on the element 
are used to obtain the two unknomis ? and n. 
equation results for 


are given, the two equations 
Eliminating n, a quadratic 


aC? + hO + c = 


0 








After is obtained from the equation above, the value for n calculated 
from the equation: 


y4 + 5(73 - ^4^ 


( 38 ) 


The resulting equation for the pressure coefficient at point i due to 
pressure at corner grid point j is: 


Cy - fj(q,np 


(39) 


6. The loads will be distributed to the structure points according to the 
location of each polyhedgon point on the structure area. The total force 
and center of force will be preserved. For each point, i, on the polyhedron 
the load factors for the structure points at locations (xj^,y 3 >, (x 2 ,y 2 >, and 
(x 3 ,y 3 ) are obtained from the determinants of matrices as shown below: 


li 


1 

A 


1 Xj yj 

1 Xo y„ 


1 Xj yj 


(40) 


21 


1 

A 


1 *1 
1 x^ y. 


1 Xj yj 


(U) 


•31 


1. 

A 


1 x^ y^ 

1 Xj yj 
1 Xj, y^ 


(42) 


7 . 


The pressures o'n the polygon points are Integrated over the area according 
to the following rules: 

a. The average pressure of all points acts over one-half the area, 
located at the center of the area. 


b. The pressure at each point acts locally over an area of A/2N. 
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, IP .^1 o' ,iy jj .p,.pi Ji 


c. The total force at the center Is divided equally among the polygon 
points. The local forces are applied at the adjacent points. 

The resulting equation for the effective polygon area coefficients is 5 


F = E A. p 
s is^^s 


(43) 


where : 




A 

s k 

2N(N- 1) 

A 

2N 

s = k 


The resulting area factor matrix is defined by the matrix product: 


[\^] = [f]'^[A][C] 


(44) 


In order to provide force vectors in three dimensions, each row of 
the matrix A is expanded to three rows by multiplication with the unit 
normal vector k, . 


2.4 GRAVITY EFFECTS 


When a steady-state acceleration such as gravity is present in a hydro- 
elastic problpm, additional terms must be added to the fundamental equations 
to account for the steady-state pressure gradient. In the fluid formulation, 
the Euler equations assume that the pressure is defined at points fixed in 
space, and the fluid particles flow across the point. In the structure for- 
mulation, a Lagrange assumption is used whereby the grid points remain 
attached to the moving system, and the forces are applied at the displaced 
location. These contradicting assumptions require formulation of additional 
matrix terms as developed below. 


2.4.1 Gravity Effects on the Structure 

A change in force on the structure is illustrated in the sketch below. 
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I** • f • 7 ™ < 


r 


/ Original Position 



Displaced Position 


-►-X 


The normal force. F^, required to support the pressure is: 

F = - A[p^ + p(g • Uj^)]n 


(45) 


The term Ap^ is included in the area matrices discussed previously. The 

second term on the right-hand side of Eq. (45) takes the form of a stiffness. 
The matrix takes the form: 


' ^ 


[u 

X 1 
1 

■ “ [K]. 

1 

F 

2 i 


l“J 


where 


(46) 


Note that the matrix is not symmetric if n^ ,< 0. This violates the fundamental 
rule that symmetric system matrices must occur for the conservation of energy. 

Another method of viewing the problem resolves the non-symmetric issue. 

If the structure moves, the total fluid weight changes as Illustrated below. 



Displaced Position 
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r.f 


The aiddUional weight on the structure, w, due to the motion is: 
w = pgAu^ = pgA(n • u) 


( 48 ) 


Since each point may move independently of the others, the increased vertical 
force must be applied locally and the force required to support the load is: 

Fz “ - P8^<" • «) * (49) 

The corresponding stiffness matrix is: ’ 



(50) 


Comparirr Eqs. (47) and (50), we observe that the lower right-hand terms are 
equal, but the off-diagonal terms are reversed. The conclusion is that each 
approach missed an off-diagonal term, and the true result is: 


[AK] « 



(51) 


These stiffness terms may be processed along with the fluid-structure area 
coefficients described in Section 2.3. The intersecting structure fluid areas 
are used to define the factor A. The displacements and resulting forces are 
assumed to be variable on the surface, dependent on the connected grid points, 
and the actual stiffness will be: 


K 


ij 




(52) 


where u^ and u^ are linear functions of the grid point displacements. These 
Integrals are evaluated in a manner similar to those developed in the area 
matrix calculations. 
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2*4.2 Free Surface Effects 


A free surface is defined as a moving boundary with no restraints. When 
gravity effects are neglected, the boundary condition, p = 0, may be enforced 
by simply applying single-point constraints (SPC) to the input which causes 
the rows and columns corresponding to zero pressure to be removed from the 
matrix equations. However, when gravity is present we must remember thit the 
pressure may not be zero since it is actually measured at a point fixed in 

space. For an upward displacement, u^, of the free surface, the pressure at 
a point defined at the surface is: 

P “ “f (53) 

(For a downward displacement, it is also convenient to use the same equation, 
measuring a fictitious negative pressure above the surface.) 

In the actual solution of the free surface points, it is convenient to 
implement Eq. (53) in the following form: 

- Ap + pgA u^ = 0 (54) 

where A is the free surface area associated with the fluid point. The terms 
in the above equation may be implemented directly into the matrix formulation. 
In effect, the free surface points are treated as though they were structure 
points, although no structural stiffness is present. The area factors A are 
Identical to the fluid /structure interface matrices defined previously in 
Section 2.3. The terms (pgA) are, in effect, positive springs providing the . 

stiffness terms, [Kfg] , for the normal displacements, uf, and causing the 
"sloshing" modes. 

Furthermore, the effects of the displacements at the free surface excite 

the fluid in the same manner as the structure displacements. The generalized, 
forces on the fluid are: 


tQf> “ - [a|]'^{u^} (5 

where [Af8] is a diagonal matrix of area factors connecting each free surface 
displacement to the corresponding pressure degree of freedom. 
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2.5 SYSTEM MATRIX SOLUTION 


The previous development has provided the basic matrix equations to 
define the fluid, the fluid structure interface, and the free surface. In 
review, these equations are: 

FLUID: 

[K^Hp} + = {0} (56) 


STRUCTURE: 

[\]{%) + IK^ + AKgHUg} - [Aj{p> = {F> (57) 

S R 

where [M ], [K®], and {F} are the conventional mi^ss, stiffness, and 
load matrices for the structure. 

FREE SURFACE 

[k|]{u^} - [Aj]{p} = {0} (58) 

2.5.1 General Formulation 

For the general case, when gravity is present, all the above matrices 
occur. The desired form of the solution matrices are: 

[M]{u> + [K]{u} = {P} ( 59 ) 

where {u} is a vector containing both structure and free surface displacements 
and {p} is the applied load vector. From Eq. (56), it is apparent that: 

{P} = - [K^]‘^[A]'’^{ii} (60) 

where [A]’’^ = [A^ j aJ] (61) 

(62) 
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and 


(u) 


U, 


Substituting into Eqs. (57), (58), and (59), we obtain; 








A 


[M] 




(63) 


M I 0 1 

s I 


oTo 


[K] 



+ AK i 

i.L„ 



(64) 


We observe that the matrices [M] and [ic] are sy mm etric, and may be processed 
as normal structure matrices. 

# 

Unfortunately, the effect of the fluid mass terms in Eq. (63) is to fill 
the mass matrir., resulting in potentially time-consuming solutions for large 
structures. However, it is typical for large structures that a reduction 

procedure is employed. Defined symbolically, this reduction may be defined 
as: 


(%} = [G]{u^> (65) 

where the vector {u^} is defined by a much smaller number of degrees of freedom 
than {ug}. Components of the vector {u^} are removed by application of con- 
att \nts through the "Guyan” reduction procedure or through a modal formulatio n 
•». the colvmms of [G] are eigenvectors of the empty structure normal modes. 
Tue structure matrices are reduced accordingly with the equations: 



[M^l = [G]^[M ][G) 

s 

(66) 


[K^'J = [G]^[K ][G] 

S 

(67) 

Equation 

(60) may be rewritten as: 



{P} = lK^j"^[Al''’{ii} 

(68) 

where 

[A]'*’ - [gV I A^] 

S * f 

(69) 


u 

a 


and 


{u} = 


(70) 




3SS 




The reduced mess and stiffness matrices are: 

! 0 


M = 


K = 


0 } 0 


+ [A][k:^][a''^] 


+ AK® ,' 
I 


where 


AK® = c’^AKgG 


(71) 


(72) 


(73) 


Note that, as the size of the matrix [A] is reduced, the evaluation of the 
matrices for Eqs. (71) and (72) will be more economical. In th. actual 
formulation, the columns of the matrix [A] may be treated as load vectors on 
the structure, and the NASTRAN reduction procedure for the load vectors may 
be applied directly. The gravity "stiffness" matrix [AK*"] may be reduced in 
the NASTRAN system with the same algorithm as the mass matrix reduction process. 

2^5.2 Non-Gravity Case 

When the effects of gravity are ignored, the free surface is constrained 
such that p = 0. In this case, the free surface pbints are removed from the 
solution vector, and the solution matrices are: 



[M] = 

[\] + [A][K^]"^[A]^ 

(74) 


[K] = 

[Ks] 

(75) 

where 

{u} = 

^-s> 

(76) 


[A]’’ = 

[Gj'^tA^]'^ 

(77) 

and [K 

j] is the matrix [K^] with 

free surface rows at»d columns removed. 


2.5.3 

Completely Enclosed Fluid 




When the fluid boundary is completely enclosed by the structure and free 
surfaces and no constraints are applied to the fluid points, the incompressible 
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fluid effects must be considered. The incompressible fluid, in effect, 
provides a constraint on the motions of the boundary such that the net flow 
into the fluid is zero. 



1 Furthermore, the fluid matrix [K^] is singular because a constant pressure 

defines zero flow. Mathematically, a unit pressure vector, defined as {l}, 
produces the result: 

j [KjHi} = {0} (78) 

Since the matri;< [K_] has a singularity of order one, a constraint must be 

i ^ 

supplied. Because of Incompressibility, we know that the total flow must be 
zero. The basic pressure-flow relationship is: 

[K]{p} = {q} (79) 

I The "average" input flow is: 

1 N 

Q = N^^i “ N ( 80 ) 

1=1 " 


where [I] is a row. vector containing unit values. 

Subtracting the average flow from each point, we obtain a flow vector, 
{Q*}, with a zero total value, where: 



{Q’} = {Q} - Q{I} 

(81) 

or 

{Q’} = [U1 ~ {I}IU]{Q} 

(82) 


= [H]{Q} 

(83) 


The pressure is obtained by removing one row and column, and solving the 
basic matrix equation in partitioned form: 






jji ^1 1.1 


' Jl. ■ I ijiM.fl lUMfWPW 


t^here a is an undetermined 
ignored, we obtain: 


constant. Noting that the {l} vector may be 
= n» 


'A - K.. QI 

J JJ 


(85) 


The coefficient a may be obtained by restricting 
have an avnrage of anro, or: fha preasnre rasult to 


P ~ N = 0 

Therefore, a is obtained from the equation: 

1 


( 86 ) 


a = - 


» (87) 

Note that Eq. (84) is completely satisfied if Eq (87) Ic « k f 
Eq. (84) and if: <!• ( 7) is substxtuted into 


N 




= 0 


( 88 ) 


and 


H 


Qi = Sq! = 0 

2=2 J 


(89) 


The first condition is satisfied 


satisfied by Eq. (82), thereby producing 


l>y (78). The second condition is 


equation. 


a unique solution to the basic 


The resultant solution pressure vector {p} is: 

{p} = 


-I iij 


ti) - 5 (mu 


(Pj) 


(90) 


Obaarva that tha prasaura transformation matrix la 5q. (90) Is Identical 
to the transpose of the flo„ trarafcrmatlon. . „hera = 


(Qj) 

%j 

(91) 

(P> 


(92) 
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and [Hj] = -^<I>Uj] (93) 

The matrix “inverse" may be written symbolically as: 

[K^]"^ = [H^]'^[Kjjl"^[H^l (94) 

Furthermore, it may be proven by examples that [K^j] may be obtained by par- 
titioning any fluid point, p^, from the matrix. If the matrix [K^] is singular 
(of order 1), the results are exactly the same regardless of the choice. 

2.5.4 Incompressible Fluid Restraint 

As described in Sectioi 2.5.3 above, the net volume change due to boundary 
movement is eliminated from the fluid inertia matrix. However, the incompressi- 
bility of the fluid requires that the volume change due to structure and free 
surface displacements be restricted. This constraint could be implemented by 
supplying a constraint equation of the form: 

AVol = T. I k..u. = 0 

i j 3^ 3 

or, in terms of the matrices: 

AVol = [I] [A]'"^{u} = 0 

For this approach, one of the displacements, , is removed from the 
matrices, redistributing its associated mass and stiffness to the other degrees 
of freedom. 

In the alternate method, we add a compressibility factor such that the net 
volume change will be small. If we define the factor, B, such that for the 
static case: 

{p} = {1} B AVol ( 97 ) 

then the forces, (f), on the structure are: 



{F} = 

[A]{p} 

(98) 

or 

{f} = 


(99) 

where 

[KJ = 

H[A]{I}HJ[A]‘^ 

(100) 
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(95) 


(96) 


•mh 


.A 


The matrix [K^] provides for the overall compressibility of the system 
when the fluid is completely enclosed. It is added to the system stiffness 
matrix and acts as a single sprin;, conn. . :ing all surface displacement degrees 
of freedom. 

The factor B may be obtained from the physical properties. The approximate 
value is: 

Pa^ 

® ~ ^ ( 101 ) 

where a is the speed of sound in the fluid. 

Extremely large values of B are to be avoided, since matrix numerical con- 
ditioning problems will result when the terms in the matrix: [K ] are orders 

c 

of magnitude larger than the structure terms. For most realistic fluid- 
structure combinations, this problem will not occur. 


3.0 TEST RESULTS 


An extensive test program was performed on the modified NASTRAN system 
during the coding effort and following the delivery of the system. The 
purpose of the program testing was to ensure correct code and validate the 
theoretical assumptions. In the first stage of check-out, problems consisted 
of simple one and two fluid element shapes in which the results could be 
hand-checked for correctness. This was followed by larger order, more 
realistic test and demonstration problems. 

The choice of test and demonstration problems had to be limxted to 
cases with known results from experimental tests and/or published analyses. 
Larger order detailed models representing the Space Shuttle External Tanks 
have also been analyzed by NASA using the program. Results of these tests 
are forthcoming from NASA. The basic test and demonstration problems 
analyzed by UAI are described below. 

3.1 HEMISPHERICAL TANK TEST PROBLEM 

3.1.1 Problem Definition 

The solution for axisymmetric sloshing and hydroelastic modes in a 
full hem spherical tank have been obtained by several methods of analyses 
[Refs. 3, 4, b and 61. Although the NASTRAN program was developed for 
general shapes, axis>-mmetric geometries such as this problem may be solved 
by modeling a wedge-shaped section with a minimum number of elements and 

grid points 

The NASTRAN model, shown in Figure 2, represents a 15“ sector, repre- 
sented by a single layer of elements. The fluid is represented by wedge, 
tetrahedron, and hexahedron elements. The structure is modeled by standard 
NASTRAN plate elements. The large sector angle and course mesh were deliver- 
ately chosen as representative of the expected modeling practices to be used 
when non-axisymmetric three-dimensional problems are generated. 

3.1.2 Comparison of Results 

Comparisons of natural frequencies for both slosh modes aad structure 
Interactions (bulge) modes are presented for the reference studies and 
NASTRAN in Table 1. Slosh mode shapes are shown in Figure 3. 






]■ 

I, 


j TABLE 1. COMPARISON OF NATURAL FREQUENCIES FOR 

} VARIOUS SOLUTIONS OF HEMISPHERICAL TANK 







Natural Frequency - Hertz 

Type 

Mode 

Ref 3 

Ref 4 

Ref 5 

3-D 

NASTRAN 


1 

0.46 

0.43 

0.43 

0.45 

Slosh 

2 

0.62 

0.62 

0.62 

0.67 


3 

0.75 

0.81 

0.76 

0.88 


4 

0.86 

1.00 

0.90 

1.13 


1 

6.69 

7.26 

6.62 

6.87 

Bulge 

2 

9.92 

11.47 

10.25 

10.99 


3 

12.59 

14.56 

12.35 

14.76 
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FIGURE 3. MODE SHAPES OF AXISYMMETRIC SLOSH MODES (1-4) 









The modal frequencies of the present analysis compare reasonably well 
with those obtained by other investigators. The higher slosh mode frequencies 
could be expected to match more closely those presented by Guyan if the num- 
ber of free surface elements was Increased from 8 (NASTRAN) to 21 (Guyan). 

This conclusion arises from the common observation that modal frequencies 
tend to decrease as the model becomes more refined by Increasing the number 
of finite elements. This explains the close match in slosh mode frequency 
for the first few modes and the divergence of results in the later modes. 

Also contributing to the differences in the NASTRAN results is the effect 
of representing the axisymmetric motions by modeling a 15“ sector whereas 
the other analyses solve the axisymmetric problem directly. The actual 
integral over the 15“ sector in NASTRAN represents a smaller area and, 
therefore, higher frequencies. This difference would be smaller for smaller 
sector angles. However, with 8 free surface elements, the present analysis 
results in gqod slosh mode frequency agreement (0-20%) . 

3.2 SRI TEST TANK 

As a further test on the performance of the 3— D analysis of a typical 
problem, a series of analyses were run on a real tank model. This actual 
model was built and tested by Southwest Research, Inc. and the experimental 
results are described in Reference 6. Other analytic results were obtained 
using the DYNAS0R axisymmetric program described in Reference 7. 

The finite element NASTRAN model is shown in Figure 4. Again, a 15“ 
sector was modeled with one layer of elements and two layers of grid points 
to solve for the axisymmetric modes. 

The mesh size was chosen such that when it was extended to a three- 
dimensional half model (12 layers) , the number of degrees of freedom (- 2900) 
would be near the maximum for reasonable running ti me . 

The effects of nearly all of the available options in the hydroelastic 
system were evaluated with the SRI model. The results are summarized in 
Table 2. The error ratios in terms of the test results are given in 
Table 3. Each of the analysis cases is described below. 

Test Results - Obtained from Reference 6. 


z 

i 


•Fluid Points \ 
/ I \ 


JjSOlO) 
(5011) 


1012 


1013 


1014 


1015 


(5012) 

\ 

Structure 

Points 

/ 

(5013) 


(5014) 


(5015) 


(5016) 


(5020) 



1021 (5021) 


1022 (5022) 


(5023) 


1029 

(5029) 


1028 (5028) 


FICUPJ: 4. FINITE ELEl-IENT MODEL OF THE SRI TANK 
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TABLE 2 


COMPARISONS OF FREQUENCIES FOR SRI TEST TANK 


Analysis Case 

Mode 

i Frequencies 

Mode 1 

Mode 2 

Mode 3 

Test Results 

495 

835 

1255 

DYNAS0R Program 

531 

807 

1179 

NASTRAN - Phase I Program 




Model A - Comp. 

519 

822 

1239 

Model B - Comp. 

516 

826 

1239 

Model B - Incomp. 

541 

828 

1240 

Model B - 1/6 Comp. 

423 

821 

1234 

NASTRAN - Phase II Program 
(Model B - Comp.) 




Direct - Not Rediced 

513 

809 

1174 

Direct - Reduced 

612 

914 

1279 

Direct - Ignore G 

539 

811 

1175 

Modal - 30 Modes 

568 

814 

1185 
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TABLE 3. COMPARISONS OF FREQUENCY ERRORS 
FOR SRI TEST TANK 


Frequency Diffsrence 
Ratios (%) 

Analysis Case — — i i~ 

Mode 1 Mode 2 Mode 3 


Test Results 

DYNAS0R Program 

NASTRAN - Phase I Program 

Model A - Comp. 

Model B - Comp. 

Model B - Incomp. 

Model B - 1/6 Comp. 

NASTRAN - Phase II Program 

Direct - Not Reduced 
Direct - Reduced 
Direct - Ignore G 
Modal - 30 Modes 


0 0 0 

7.3 -3.35 -6.1 

4.85 -1.56 -1.28 

4.25 -1.08 -1.28 

9.3 -0-83 -1.20 

-12.5 -1.68 -1.67 

3.7 -3.1 -6.5 

23.6 9.5 

8.9 -2.9 

14.8 -2.5 


-6.4 

-5.6 









DYNAS0R Program - The tank was modeled and run at MSFC on the DYNAS0R 
Program. 

NASTRAN - Phase I Program - The first system delivery contained limited 
options and a crude method of calculating area coefficients. Nc over- 
i- pplns structure/fluid elements were allowed. All runs were made 
using the direct formulation method with no matrix condensation. 

Model A ~ Compressible — This model was generated by simply con— 
^®^bing the DYNAS0R data to the NASTRAN format. The mesh was similar 
to that shown in Figure 4 except that only fotir— sided elements were 
used. The compressibility factor was obtained from the properties 
of water. 

2) Model B - Compressible - This was the basic test case using the 
model shown in Figure with overall compressibility of water. The 
second and third modes were excellent but the first mode was 
suspiciously high. 

3) Model B - Incompressible - The incompressible option was used in 
this model to determine its effect. The first mode became worse but 
the second and third modes were hardly affected. 

4) Model B - 1/6 Compressibility - The compressibility factor was 
divided by a factor of 6. The first mode frequency became lower 
than the test results with no change in the second and third modes. 
Tills indicated that fluid compressibility had affected the test 
results. The displacement! , in the first mode were primarily bulging 
of the tank, with nearly uniform vertical motion of the free surface, 
resulting in a net total pressure over the fluid interface. The 
second and third modes contained little net free surface motion and 
net pressure was small. The actual compressibility of the water in 
the test was probably lower than the theoretical factor due to 
aeration of the water during the vibration testing. 

NASTRAN - Phase II Program - The final delivered program contained more 

accurate area factor calculations and the complete set of user options. 

The tests given below were run on this version. For comparison wich the 
P^^®lin>lnary version. Model B with the calculated compressibility was used 
as the basic model. 
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Direct Not Reduced - The direct method without matrix condensation 
was used to compare results with the Phase 1 program. Results for 
the first mode were improved but the second and third modes became 
slightly worse. 

6) D irect - Reduced - In this case the solution matrices were reduced 
from ~257 degrees of freedom to ~60 degrees of freedom to represent 
only shell displacements at every other point. This reduction would 
be equivalent to reducing the three-dimensional model to ~300 degrees 
of freedom for eigenvalue extraction. All modes incre->sed in frequency. 

Direct - Ignore Gravity - The gravity effects were removed from the 
problem which reduced the solution size and the running time. The 
small change in results indicates that this is an efficient method for 
obtaining structixre interaction modes. Low frequency slosh modes may 
not be calculated with this method. 

8) Modal — 30 Modes - The modal formulation was used in this problem 

to reduce the structure matrices to 30 modal coordinates representing 
the modes of the empty structure. The error in the first hydro- 
elastic mode was due to the fact that its shape was not well repre- 
sented by the mode shapes of the empty structure. Only three of the 
30 empty structure modes participated to any extent in the first mode 
of the combined fluid and structure systems. 

3.3 TEST RESULT COMMENTS 

From the experience of running the test and demonstration problems, 
several conclusions may be made regarding the NASTRAN hydroelastic system. 

These are listed below. 

1. Accuracy of the system was better than expected for the mesh sizes 
used in the demonstration problems. With only linear elements and 
averaged area factors representing the fluid, three good slosh modes 
were obtained from only eight degrees of freedom. It appears that 
the accuracy for hydroelastic modes is limited more by the existing , 
NASTRAN structure elements than by the fluid formulation. Results 
indicate that 15“ sectors are adequate for a cylindrical or spherical 
shaped fluid model. 


1 
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The results were relatively insensitive to modeling procedures. On 
each of the problems, different methods of subdividing the fluid 
space into elements were tested. For similar mesh sizes, the changes 
in results were insignificant. 

The use of either Modal Formulation or Guyan reduction to condense 
the structural degrees of freedom tends to increase the natural 
frequencies of the system. For the relatively small demonstration 
problems, their effects on execution cost were small. However, the 
hydroelastic formulation produces dense solution matrices. Large 
order problems will require one of these reduction methods. 

The drawback to the Guyan reduction method is that it poorly 
represents the motions of curved surfaces with uniform loads. For 
best use of this method, all displacements normal to the surface 
should be retained. In-plane displacements may be omitted without 
affecting results since they are not connected to the fluid mass. 

Modal formulation is best used when the empty— structure modes 
are similar to the coupled fluid modes. Some of the low frequency 
combined— system modes do not occur as the lowest modes for the empty 
structure. Thus, for some cases, many modal degrees of freedom may 
be required to produce accurate results. 

Although free-surface gravity effects are necessary to obtain pure 
sloshing modes, their effect on the hydroelastic modes is small for 
most problems. The alternate method of constraining the free surface 
pressures to zero is more efficient and requires less data input. 

The overall compressibility factor used in the new method provides 
a simple, efficient manner of treating enclosed fluids. It will 
produce more accurate results for very stiff tanks such as those 
used in the SPI demonstration problem. 
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